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Studies of fixation dynamics in Markov processes predominantly focus on the mean time to absorp¬ 
tion. This may be inadequate if the distribution is broad and skewed. We compute the distribution 
of fixation times in one-step birth-death processes with two absorbing states. These are expressed in 
terms of the spectrum of the process, and we provide different representations as forward-only pro¬ 
cesses in eigenspace. These allow efficient sampling of fixation time distributions. As an application 
we study evolutionary game dynamics, where invading mutants can reach fixation or go extinct. We 
also highlight the median fixation time as a possible analog of mixing times in systems with small 
mutation rates and no absorbing states, whereas the mean fixation time has no such interpretation. 
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I. INTRODUCTION 

If a group of mutants is introduced into a popula¬ 
tion of wild-type individuals, how long does it take for 
the population to reach a homogeneous state? This 
is a typical question asked in population genetics and 
evolutionary biology. It can be addressed using the 
theory of stochastic processes and approaches from 
statistical physics. 

In its simplest form, evolutionary dynamics is mod¬ 
eled as a Markov process in a population of N indi¬ 
viduals of two types [1]. Fixation occurs when the 
process arrives at one of its absorbing states. As the 
time to fixation is itself a random variable, only the 
computation of the distribution of fixation times pro¬ 
vides a complete answer to our opening question. The 
majority of existing studies avoid this mathematical 
challenge and instead focus on calculating only the 
mean fixation time mu. While the first moment can 
provide a good indication of the outcome in some cir¬ 
cumstances, this approach can be insufficient when 
the distribution of fixation times is broad [SI IS]. As 
we show in this paper, such scenarios arise in exam¬ 
ples from evolutionary game theory. 

Although the master equation describing the birth- 
death dynamics is linear, calculating fixation time 
distributions is more intricate than one may initially 
think. Nested expressions for all moments of fixation 
times are known Pl] and from these the distribu¬ 
tion can in principle be constructed recursively up to 
arbitrary precision. However, this approach does not 
provide a simple closed-form solution or a means of 
efficiently sampling from the arrival time distribution. 
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An alternative approach is to diagonalize the linear 
operator of the master equation and to carry out the 
analysis in eigenspace. This leads to a theorem at¬ 
tributed to Karlin and McGregor PHo], which states 
that arrival times can be written as the sum of inde¬ 
pendent, exponentially distributed random variables 
parametrized by the eigenvalues of the master equa¬ 
tion. Alternatively stated, the arrival time distribu¬ 
tion is given by the convolution of the exponential dis¬ 
tributions from which the independent random num¬ 
bers are sampled, i.e., it is a phase-type distribution 
m- This theorem has been discussed in numerous 
sources in the probability theory literature nni nu¬ 
lls]. The discussion of these matters is usually very 
terse, and not easily accessible to physicists or re¬ 
searchers in adjacent disciplines. Researchers in the 
theoretical biosciences are only recently beginning to 
use these ideas for the purpose of model reduction 
OEI]. Existing results are limited to specific ini¬ 
tial conditions and types of birth-death chains, and 
a clear understanding of the analysis in eigenspace is 
lacking. 

We consider one-step birth-death processes with 
two absorbing states and a general initial condition. 
This model describes the invasion (or extinction) of 
a number of mutants in a wild-type population. It is 
a generalization of existing studies, which focus on a 
single absorbing boundary and a specific initial condi¬ 
tion |9| . The purpose of our work is to provide a sys¬ 
tematic and explicit closed-form solution for fixation 
time distributions, and a physical interpretation of 
different representations in eigenspace [HIITIIIH]. By 
generalizing the Karlin-McGregor result to account 
for arbitrary initial conditions and multiple absorb¬ 
ing states, we show that these different representa¬ 
tions can be traced back to one common origin. To 
generate samples from the fixation time distribution 
it is sufficient to simulate forward-only processes in 
eigenspace. This provides effective model-reduction 
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one-step birth-death process in real space 
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FIG. 1. (Color online) One-step birth-death process in a 
population of N individuals. The variable i denotes the 
number of invading mutants. The states 1 = 0 (extinc¬ 
tion) and i = N (fixation) are absorbing. Birth rates are 
labeled bi and death rates di. 


schemes. 

We use our results to relate fixation processes to 
the equilibration dynamics of evolutionary systems 
with mutation (and hence with no absorbing states). 
These equilibration processes are often characterized 
by the so-called mixing time, the time it takes the 
system to come within a set distance of its stationary 
distribution [53] . We thus discuss an appropriate 
analog of the fixation time in the limit of small muta¬ 
tion rates, and outline the limitations of this relation. 

The paper is organized as follows: In Sec.jT^we de¬ 
scribe the model used and the principal idea behind 
the work. In Sec. m we outline the mathematical 
procedure that takes us from the model to the arrival 
time distributions. Details of the calculations can be 
found in the Appendixes. In Sec. jlVj we apply our 
results to specific examples of evolutionary games, 
confirming both the accuracy and efficiency of our 
method. In Sec. we discuss the relationship be¬ 
tween the fixation time and mixing time in the limit 
of small mutation rates, before drawing conclusions 
in Sec. |Vl| 


II. MODEL DEFINITION AND MAIN IDEA 

We study a one-step birth-death process with 
states i = 0,1,... ,N and characterized by the birth 
and death rates bi,di {i = 1,..., N — 1), as illustrated 
in Fig. [2 This describes a population of constant size 
N with i individuals of the mutant type and N — i 
of the resident wild type [Tj. The states i = 0 and 
i = N are absorbing in the absence of mutation, and 
so the dynamics will end at one of these two states 
eventually. Our objective is to calculate the distribu¬ 
tion of arrival times at these absorbing states for a 
general starting point *o- 

The outcome of our analysis (detailed below) are 
the reduced forward-only processes in eigenspace il¬ 
lustrated in Fig. which have the same arrival time 
distribution as the original birth-death process. The 
forward jump rates are determined by the eigenvalues 
Aq, of the original birth-death process. The schematic 
in Fig. Ha) shows multiple forward-only channels, 
where dashed arrows indicate that some steps are to 
be skipped. Figure [^b) shows a single forward-only 



FIG. 2. (Golor online) Different representations of the 
arrival process in terms of forward-only processes in 
eigenspace. The \a are (absolute) eigenvalues of the pro¬ 
cess in Fig. each arrow represents an exponential pro¬ 
cess with the rate indicated, (a) shows a set of alter¬ 
native reaction channels. In each run one channel is cho¬ 
sen with appropriate probability. Transitions indicated by 
dashed arrows are skipped (zero time), (b) shows a single 
forward-only chain, in which the final state can be reached 
directly from some of the intermediate states. Both rep¬ 
resentations are equivalent and generate samples from the 
arrival time distribution of the process in Fig.[2 The case 
shown here is for arrival at N, starting from io = 3 in the 
original space. 


channel. Long arrows indicate direct jumps to the fi¬ 
nal eigenstate. In both cases, the number of possible 
paths is determined by the initial condition and the 
final state (f = 0 or i = N) of the original birth-death 
process. 

Arrival time samples of the original process are gen¬ 
erated from Fig. Ha) in the following way: One of the 
channels is chosen with probability determined by the 
birth-death rates, the initial condition, and the ar¬ 
rival state of the original real-space process. After a 
channel has been selected, the clock is started and the 
forward-only process of the channel is traversed. The 
clock is stopped when the final state in the schematic 
is reached (absorption). 

Arrival time samples are generated from Fig. [^b) 
by traversing the forward-only chain. The quanti¬ 
ties 1 — Fa are the probabilities to reach absorption 
directly from the intermediate eigenstate a. These 
probabilities are again determined by the birth-death 
rates, the initial condition and the final state of the 
original process, and explicit formulas can be found 
in Appendix |C|2| 


III. ANALYSIS 

To derive these results it is convenient to focus on 
the interior states f = l,...,A — lin the birth-death 
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process in Fig. The probability to be in state i at 
time t, Pi{t), satisfies the equation p = A • p with 
formal solution p(t) = exp(At) • p(0). The matrix A 
describes the transient states and has elements = 
— {bi+di), = di+i and ai^i-i = 6i_i. The initial 
condition is Pi(0) = < N — 1). Probabil¬ 

ity continuously flows into the two absorbing states. 
Hence all eigenvalues of A are negative, and so we fo¬ 
cus on —A and its eigenvalues. The matrix exponen¬ 
tial can be evaluated by Laplace transformation, and 
subsequent back transformation allows one to calcu¬ 
late PoKo(^) = dipi{t) and pN\io{t) = bN-iPN-i{t). 
Mathematical details can be found in Appendix 
Up to normalization these are the conditional arrival 
time distributions at states 0 and N respectively. For 
the remainder of the paper we focus only on absorp¬ 
tion at state N. Analogous expressions for absorp¬ 
tion at 0 can be found in the Appendixes. We find 
the following expression for the probability flux (per 
unit time) into state N 

PNHo{t) = ^^E^_,.R,,_u ( 1 ) 

where Bi^ = ^ ~ det(—A) and V'io is the 

determinant of the (*o — 1) x (ip — 1) top-left sub¬ 
matrix of —A (denoted as —We have intro¬ 
duced Ei = ■ ■ ■ * where the i(f) are 

exponential distributions with the eigenvalues of —A, 
Aq, > 0, as parameters. The symbol * represents a 
convolution. The object is of the form 

Ri = (^ + * ■ ■ • * {S + ^6'^ , (2) 


where ?/a > 0 are the eigenvalues of the sub-matrix 
_A(*o-i)_ jn these expressions 6 is the usual Dirac 
distribution, and S' is its derivative, as described 
in Appendix |A|1[ We can identify the prefactor 
BioPio/R = 4’N\io the probability that the origi¬ 
nal system gets absorbed in state N (as opposed to 
0 ). 

The effective dynamics shown in Fig. [2| are obtained 
by evaluating the convolutions in Eq. (111). To illus¬ 
trate this it is useful to first consider the convolution 
of the exponential distribution £^^'>{t) with an ob ject 
of the form d{t) + (y > 0). In Appendix A 

we show that 


£P>*(S + y-^S') 


%)+ 1-^ £(^)(t) 


■ (3) 


Assuming A/y < 1 (which will be the case throughout 
our analysis) this describes a convex combination of 
a point mass at zero and an exponential distribution. 
For samples, set t = 0 with probability A/y, otherwise 
draw t from £^^^t). 

To arrive at the dynamics depicted in Fig. [^a), 
each of the Iq — 1 terms of the form A -|- in 

Eq. Q is paired up with a separate exponential, as 
described in Appendix |C|1| This creates a total of 


2 * 0-1 possible forward-only channels with up to *o ~ 1 
exponential steps skipped, as shown in Fig. [^a). To 
arrive at the dynamics shown in Fig. ib), the to — 1 
objects of the form 5 + y^^S' are successively con¬ 
voluted with the full chain E^-i from the right, as 
described in Appendix |C|2| This leads to io chan¬ 
nels, in which 0,1,..., Jq — 1 exponential steps are 
skipped. This can be illustrated as a single forward- 
only channel, with appropriate rates of jumping from 
intermediate eigenstates to absorption. By consider¬ 
ing these forward-only processes we arrive at closed- 
form expressions for the arrival time distributions. 
For arrival at state N, the distribution is given by 
(see Appendix |C|3[) 


N-l 


PN\ioP') — Pio ^ ^ 




io — 1 

n iVl ~ ^a) 

7=1 

N-l ' 

n (A/3 - Aa) 
/3=1 




(4) 


The representation shown in Fig. [^a) corresponds 
to the picture obtained for a restricted set of processes 
by probabilistic methods in Ref. mi. On the other 
hand. Fig. [^b) reflects the findings of Refs. mnn], 
derived from the construction of intertwining pro¬ 
cesses. Our analysis shows that these different de¬ 
compositions originate from one common structure, 
Eq. Q. The explicit schemes in Fig. [^provide a 
computational method to generate samples from the 
arrival time distribution efficiently, for example by 
carrying out simulations of these forward processes 
using the Gillespie algorithm [24]. It is important to 
keep in mind that the eigenstates shown in Figs.j^a) 
and Hb) cannot be mapped one to one to the states 
in Fig. The equivalence of the real and eigenspace 
representations only holds on the level of arrival time 
statistics. 


IV. EVOLUTIONARY GAMES 

As an application of this theory we now consider 
examples of evolutionary dynamics with frequency- 
dependent selection [U |25l [26| . Such models are used 
to describe the interaction of invading mutants (A) 
and resident wild-type individuals (B). These sce¬ 
narios can be formulated as evolutionary games [T|. 
We focus on a 2 X 2 normal form game. 



A 

B 


_ 2—1 p 1 N—i c 


A 

R 

S 



(5) 

B 

T 

P 


i rp 1 N—i—\ -p 
~ N-1^ N-l 



where 7r^(t) and 71 ^( 1 ) are the expected payoffs in 
a population of i individuals of type A and N — i 
individuals of type B. For this example we as¬ 
sume a pairwise comparison process, leading to birth 
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FIG. 3. (Color online) The conditional fixation time dis- 
tribntions at i = A^ for different evolntionary games. Main 
panel: Lines show results from the theory [Eq. Q], sym¬ 
bols are from simulations (10® runs per gam^. The 
mean fixation time (arrows) is not a good description 
of the distribution. Inset: Open markers/dashed lines 
show compnter time needed to obtain arrival time dis- 
tribntions from simulations, full markers/solid lines are 
for the semi-analytical approach and indicate the poly¬ 
nomial scaling. Direct simulations can only approximate 
the exact arrival time distribntion to a given accnracy 
(see Appendix C|4 for details). (A = 100, io = 10, /3 = 0.1. 
Coexistence game: R = P = 1.0, S — T = 1.5; Coordina¬ 
tion game: R — P=1.5, S = T =1.0; prisoner’s dilemma: 
R = -S = 0.5, r = 1.0, P = 0.0). 


and death rates given by bi = g[+ATr{i)]i{N — i)/N 
and di = g[—ATT{i)]i{N — i)/N, respectively, where 
A7r(i) = 7r^(i) — TrB{i) and g{z) = (1 -h [3z) /2. The 
parameter /? > 0 is the so-called intensity of selection 
[DimiiH]- The parameters R, S, T, P specify the in¬ 
teraction, and we consider three possible scenarios: 
coexistence {T > R and S > P, stable heterogeneous 
population); coordination [R > T and P > S, un¬ 
stable heterogeneous population); and the prisoner’s 
dilemma {T > R> P > S, stable homogeneous popu¬ 
lation) [1]. 

As shown in Fig. arrival time distributions can 
be broad and skewed, such that the mean fixation 
time contains only limited information. Figure also 
demonstrates the computational benefits of our for¬ 
malism. Evaluating the distribution in Eq. ([^ re¬ 
quires O (N^) operations as the spectra of —A and 
_Abo-i) 

are needed. Generating arrival time distri¬ 
butions from simulation of the original birth-death 
process instead can take exponentially long. Eor the 
coordination game and the prisoner’s dilemma, fixa¬ 
tion at N is rare and the bottleneck of direct simula¬ 
tions is the limited sampling. Eor coexistence games, 
fixation times are exponentially long in N. This im¬ 
pedes accurate simulations. Direct numerical inte¬ 
gration of the master equation would suffer from the 
same problem. 


V. EQUILIBRATION PROCESSES IN 
SYSTEMS WITH MUTATION 


We now consider birth-death processes without ab¬ 
sorbing states by adding mutation occurring at a rate 
u <C 1, such that bo = O (u) and = O (u). All 
other transition rates depicted in Fig. [ll are O (u°) 
and are only affected at sub-leading order by u (ex¬ 
act transition rates can be found in Appendix |D|2[ ). 
The timescale of the dynamics is characterized by 
the so-called ‘mixing time’, tmix- This is the time 
taken for the probability distribution, P(t), to come 
within a specified distance of the stationary dis¬ 
tribution P®*, i.e., finix is the first time at which 
d[P(tniix)j P*^*] = £■ The distance between distri¬ 
butions P and Q commonly used in this context is 
d(P, Q) = jPi - Qi|/2 with e = 1/2 [m |23]. 
Using our results we can determine if and when there 
is a correspondence between the mixing time and the 
fixation time. 

For very small mutation rates (0 < uN <C 1) the 
stationary distribution is of the form Pf* « (1 — 
o')<5i,o + where a = O (u°) (see Appendix D 2). 

In the strict absence of mutation {u = 0) the sys¬ 
tem reaches fixation, so its terminal distribution, $, 
is of the form = (1 - 4>N\io)Si,o + (l>N\io^i,N- It is 
clear that these two distributions are different; P®‘ in 
systems with rt > 0 is independent of the initial con¬ 
dition. Thus there is no obvious connection between 
fixation times and mixing times in the limit u —>■ 0. 

However, equilibration in many systems with rare 
mutations is a two-step process; the system first 
reaches an intermediate distribution that is depen¬ 
dent on the initial condition, before ‘leaking’ on a 
longer timescale into the final stationary state [53]. 
This is analogous to the quasi-stationary distribution 
before reaching absorption in systems without muta¬ 
tion [T3] . Our analysis suggests that this intermediate 
distribution of systems with 0 < uN 1 is close to 
the terminal distribution of the system with u = 0, 
and that both systems initially evolve along similar 
trajectories (see Appendix |D|2[ ) . This can be seen in 
Fig.j^where we represent the probability distribution 
as single points in the (Pq, PAr)-plane. The distribu¬ 
tion first approaches the terminal distribution $ be¬ 
fore slowly converging to the stationary distribution. 
The most appropriate analog of fixation times in sys¬ 
tems with small mutation rates is thus the time to 
reach this intermediate distribution, not the time to 
stationarity (mixing time). 

The mixing time can related to the fixation time in 
a different way: The probability distribution initially 
satisfies (i[P(t),P®‘] = d[P(t),4>], which is a conse¬ 
quence of the distance measure used (solid lines in 
Fig.lU). This equivalence hol ds w hile Po(t) < 1 — cr 
and PAr(t) < cr (see Appendix |D|2[ ). Prior to the first 
time that this condition is violated, P(t) is approxi¬ 
mately the same as in the system without mutation 
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FIG. 4. (Color online) Approach to the stationary dis¬ 
tribution for a system with small mutation rates. Reac¬ 
tion rates are given in Eq. ( |D6[ ). Dots show the trajec¬ 
tory of [Po“' (t ), (t)] (found from numerical integra¬ 

tion of the master equation with mutation). The prob¬ 
ability quickly approaches the fixation distribution, 
before slowly leaking to the stationary state, For 

0 < io < the trajectory starts at (0,0) and leaves the 
shaded area at time t*. For any point inside the shaded 
area, the distance to the points (1 — <(>jv|io,'('iviio) 

(1 — (7, a) in our metric (solid lines) are equal. Payoff pa¬ 
rameters are R = P = 1.25 and S = T = 0.75, which 
corresponds to a coordination game. Remaining parame¬ 
ters are P = 0.1, N = 100, io = 35, and u = 10“®. 


where we have Pr(lflx < t) = 1 — (i[P(t),#] for the 
cumulative fixation time distribution. The condition 
d[P(<), $] = 1/2 therefore translates into the time at 
which half of the samples have reached fixation, and 
provided the above conditions hold there is a corre¬ 
spondence between the mixing time and the median 
fixation time. This is illustrated in Fig. where we 
consider the example of a coordination game with 
a symmetric stationary distribution [53]. If P®* is 
asymmetric and a (1 — cr) (or vice versa), then the 
condition P/v(t) < cr is violated after a short period 
of time. In such cases the above correspondence may 
not hold for the median fixation time, but only for 
a lower percentile. This is the case for the example 
shown in Fig. 



Arrival time PDF Mutation rate, u 


FIG. 5. (Color online) Correspondence of mixing time 
and median fixation time for small mutation rates, (a) 
shows the unconditional fixation time distribution for the 
coordination game (u = 0). (b) depicts the mixing time 
for u > 0 and e = 1/2 (from numerical integration of the 
master equation [23)1. Reaction rates and parameters are 
as in Fig. [^ 


rival time distributions reduces to simulating these 
forward-only processes, or equivalently evaluating a 
finite sum of exponential random variables, turning 
our results into an effective tool for model reduction. 
The compact structure of the forward-only processes 
allows us to derive exact, closed-form expressions for 
the arrival time distributions of the original process 
in terms of its spectrum. As we have demonstrated, 
the numerical evaluation of these expressions is an ef¬ 
ficient polynomial-time method to obtain full arrival 
time statistics. We have also established a link be¬ 
tween equilibration times in systems with small muta¬ 
tion rates and the median fixation time in absence of 
mutation in some circumstances. Our work advances 
the understanding of the mathematical structure un¬ 
derpinning the dynamics of fixation. We have placed 
existing representations for simpler cases into a wider 
and more coherent context |141 flTl [18] . Neverthe¬ 
less, there are fundamental open questions. Claims 
of probabilistic interpretations of Karlin and McGre¬ 
gor’s theorem have been made [m dn, but in our 
view this picture is still incomplete. We would argue 
that a full probabilistic interpretation of the represen¬ 
tations in eigenspace is only reached when each time 
step of the forward-only process can be constructed 
directly and uniquely from realizations of the origi¬ 
nal process alone. Whether or not this is possible is 
unclear. 


VI. CONCLUSION 

In summary, we have constructed eigenspace repre¬ 
sentations that capture the full arrival time statistics 
of one-step birth-death processes. The mapping into 
eigenspace has a clear interpretation as forward-only 
exponential processes. Sampling of the original ar- 


Appendix A: Mathematical background material 

1. Dirac distribution and its derivative 

The Dirac-(5 distribution has support {0}. It 
can be written in its Fourier representation, S(t) = 
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FIG. 6. (Color online) This figure is as described in Fig.[^ 
except the dynamics now follow the prisoner’s dilemma 
scenario. At time t* we have d{t*) « 0.65, which means 
the mixing time and the median hxation time do not cor¬ 
respond in this example. Parameters are R = —S = 0.5, 
T = 1.0, P = 0.0, /3 = 0.025, N = 100, io = 50, and 
u = 10 “®. 


The Laplace transform is obtained as follows 




dt 


10 - 


= A 


-(s+A)t 


(A4) 


This integral is only convergent in the region Re(s) > 
—A. Within this region we have 

A 


£ 




5 “h A 


(A5) 


3. Laplace transform of an object S{t) + z 

We now show that the Laplace transform of S{t) -h 
z~^d'{t) is 1 -I- z~^s. The object 6'{t) is the der iva- 
tive of the Dirac-J distribution 6{t) (see Appendix |A|l| 
above) and z > 0 is a constant. We have 

POO 

C[6{t )= / [S{t) dt 

Jo~ 

poo 

= / e-^^S{t)dt^z-^ 

Jo- 

poo 

+ z~^s / e“®*(5(t)dt 

Jo- 

= l + z-'^s. (A6) 


duje^'^*. The distributional derivative, 5' can be 
conveniently defined by its Fourier transform as well 
P5] . It has the form 

/ OO 

(zw)e*‘"‘ dw. (Al) 

-OO 

When this is convoluted with a test function, /(t), 
with infinite support one obtains (after integration 
by parts) 

/ OO 

5'{t-T)f{T)dT = f[t). (A2) 

-oo 


then one hnds 


poo 

/ 5'{t- T)g{T)dT = g\t) + g(0)5{t). (A3) 

Jo 


This expression has no singularities, and thus the re¬ 
gion of convergence in terms of s is the entire complex 
plane. Note that we explicitly define the lower inte¬ 
gration limit as 0“ to include the 5 function in the 
integral, and we have used lim(_^o- = 0- 


4. Convolution of an exponential distribution 
with an object 5{t) -\- 

The convolution of an exponential distribution 
with an object of the form S{t) + z~^d'{t) (for 


(A2) 

constant z), as shown in Eq. is 




>0, 

pOO 

= / [S(t — t) + z~^d'{t — 

Jo 

t)] dr 


= Ae-^‘ -b \z-^Sit) - A^z-^e-^* 


(A3) 


(A7) 


2. Laplace transform of an exponential 
distribution 

We frequently use the Laplace transform of an ex¬ 
ponential distribution in our subsequent derivation. 
This is a standard result, but it is useful to include it 
here. We consider and exponential distribution with 
parameter A > 0, such that = Ae“^‘ (i > 0). 


Appendix B: Calculation of absorption time 
distributions via Laplace transforms 

1. Laplace representation 

As mentioned in the main text it is convenient to 
focus on the states i = 1,...,A^—lof the birth-death 
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process shown in Fig. 1 of the main paper. The 
dynamics of these states is given by p = A • p, where 
A is an (A^—1) x (A^—1) matrix, and where Pi{t) is the 
probability that the system is in state i at time t. We 
note that this is not a probability-conserving master 
equation, as probability mass continuously leaks into 
the absorbing states. We will use the notation pi 
(lower case) when we discuss the restricted system 
(with Pi(t) < 1), and we will write Pi(t) (upper 

case) when we discuss the full system, i = 0,..., A^. 
For the latter one always has (^) = 1 at all 

times. 

The formal solution of the equation p = A • p, 
restricted to sites i = l,...,A^ — 1, reads 

p{t) = exp(At) • p(0). (Bl) 

We can take the Laplace transform and write the 
matrix exponential in the resolvent form 

p(s) = (si-A)"i •p(0). (B2) 

We have here written p(s) for the Laplace transform 


£[p(t)]. We consider initial conditions of the form 
Pi(0) = Si^ig (1 < iQ < A^ — !)■ Our strategy is to 
compute Pi{t) and PN-i{t)^ and from these the rates 
dipi{t) and bN-iPN-i{t), with which probability ar¬ 
rives in the absorbing states, can be obtained. Thus 
we are interested in 


PiKo(s) = [(sl-^) 

PN-l\^o{s) = [(sl-A)-1]^_^_.^. (B3) 


The (i, j)th element of the inverse of an invertible 
matrix B is given by [B“^]y = Cj^i/ detB, where 
is the (j, i)th co-factor of B. Thus we can write 


Pl|io(s) = [(sl-^) \ 


and likewise for the {N — 


1 

det(sl —A) 
1, zo)th element. 


a- 




(B4) 


2. Calculation of co-factors 

The CO- factor Cig^i is found from dropping column 
1 and row zq from si — A and evaluating 



-d2 0 

S -|- 02 —do 

0 








—62 s - 1 - 03 

—d4 

0 






ao,i = (- 1 )*"+' 

0 

-bio-2 

s + dio-l 

dig 

0 



, (B5) 



0 

0 

-big 

S + Oig + l 

-dio+2 

0 






0 

0 

— bN-2 

s + aTv —1 



where ai = bi + di is introduced for compactness. Using Laplace’s formula, this can be written as 


Q„.i = (-i)*“+M];i-d. 


Vi=2 


S + Clig + l —dig+2 0 

~bio + l S -|- aig+2 —dig+3 0 


= det (si-A(7 v-*o-i)) 

io N-ig-l 
— di (s -f Xq,). 


i—2 a—1 


0 —&V-3 S-I-aAr_2 —dv-i 

0 —&V-2 s + aju-i 


(B6) 


The matrix A(jv-io-i) consists of the rows and 
columns io -f 1, • ■ •, A^ — 1 of the matrix A, i.e., it 
is the bottom right (A^ — io — 1) x (A^ — io — 1) sub¬ 


matrix of A. The matrix —has eigenvalues 
Xa > 0 {a = 1,..., N — io — 1) and determinant 
det (-A(7v-io-i)) = Xio = Xc- 

For the (io, A^ — l)th cofactor we have 



0 


/N-2 


= (-i)*o+^-i Yi -h 


s + oi —d 2 0 

—hi s + a 2 —d^ 


0 htQ — S ^ “t“ ^zq —2 di^ — i 

0 -bio-2 s + Gio-I 


rN-2 


= ( n ^ 

\i—io / 

N-2 20-1 

= n n ('5+^“)- 


(io-l)' 


2=20 Oi — l 

The matrix consists of rows and columns 

1,... ,io — 1 of the matrix A, i.e., it is the top-left 
(io — 1) X (ig — 1) sub-matrix of A. The matrix 
_A(*o-i) lias eigenvalues yo, > 0 (a = 1 ,..., io — 1) 
and determinant det (—= tpig = 2/«- 

3. Arrival time distributions 

Putting things together we have 

2o N — io — 1 N—1 


5 = 


£•0 -IV —ty —j. 2V — j. ^ 

]Jdz + n s +A/?’ 


i—2 Q;=l 

N-2 20-1 


N-1 


1 


PAr-l|zo(s) = n n + n ^ThT' 

z —zq q—1 /3—1 ^ ^ 


Here we have used det(sl — A) = 11^=1^ (s + '^/s); 
where A /3 > 0 are the eigenvalues of —A. Using 


(B7) 


Po{s) = dipi{s) and Pn{s) = bN-iPN-iis), we ob¬ 
tain the Laplace transforms of the absorption time 
distributions at sites i = 0 and i = N, respectively. 
They are given by 


io M-zo-l 


N-1 


■Poizo(s)=n (s+a;a)n 

2—1 a—1 0—1 

^ N-1 20-1 N-1 


S + A/3 ’ 


1 


-paiizo(s)= n n n ^TAT- 

i—io a—1 /3—1 ^ ^ 


4. Time representation of solutions 


Combining Eqs. (A51, (A 61 , and (B91, we can write 


Zo Af—Zo —1 


N-1 


^oiz„(s)=n'^* n ^[xaSit)+6'{t)]i[ 


i—1 a=l 

N-1 20-1 


A/3 


PN\iois) = b, C[ya5{t)+5'(t)\ J| 

2 = 2 o CX. — 1 0 — 1 


A/3 


(BIO) 


Using t he fa ct that C ^ [/(s) • g{s)] = f * g, we can perform the inverse Laplace transform of the expressions 
in Eq. plOl). We find 


-Po|zo(^) — 


-PAr|zo(^) = 


20 N-io-1 

n n 

2 = 1 Q; = l 

* N-1 

n 

p=i 

N-1 20-1 

n n 2 /o 

i—io Q!=l 

N-1 

n ^/3 


fAl) * ... * gi^N-l) * ^,5 _|_ 2;^ * . . . * _|_ X0_-^_N) J 




(Bll) 
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This is the expression shown in Eq. Q. We iden¬ 
tify the prefactors in square brackets as the fixation 
probabilities 4>o\io 4>N\io^ respectively. 

5. Symmetry of the fixation time distributions 


tional arrival time distribution at state i = N given 
io = 1. This symmetry has been known for the mean 
fixation time HEO], and it was recently shown that 
the correspondence holds for the full distribution m- 
Our approach offers an alternative way to obtain this 
intriguing result. 


By choosing the initial conditions *o = ~ 1 and 

zo = 1, the expressions in Eq. (Bill can be reduced 
to 


-Po|Af-i(^) — 


N-l 

n 

. 2=1 




= 4>o\n-iEn-i, 


PN\l{t) = 


'N-l 

n 




= (B12) 

where Eg = £Oi) * ... * From this we conclude 

^0|iV-l(^) _ PNllit) 


</'0|Af-l 4’N\1 


that is the conditional arrival time distribution at 
state i = 0 given zq = fV — 1 is equal to the condi- 

I 


Appendix C: Computing eigenspace 
representations 

As the convolution operator is commutative, we 
can order the convolutions in Eq. (Bill in any way. 


The exponential distributions and the objects of the 
form S{t) + in Eq. (Bll) can hence be com¬ 

bined in multiple ways. 


1. Convolutions I: Pairing 5 + 5' with 

individual exponential distributions 


(B13) In this section we choose to couple with 

5{t) -\- for the purposes of Po|io(^)> and with 


5{t) + yaS'{t) when we calculate Carrying 


out these convolutions in Eq. (Bll), we arrive at 


+ ^ S +1^1_£-(Aig+l) 

XjV — iQ — 1 \ — iQ — 1 / 

^Af|io(^) = ^N\io X * ■ ■■ * 

_ Uio — l \ Vio — l ) 


Xi \ Xi J J 

^ (Cl) 

. yi \ yi J 


We stress that the objects 6{t) + x~^S'{t) [or S{t) -I- 
z/“^5'(t)] can be paired with any of the exponential 
distributions. We chose to match these at the end of 
the exponential chain so that the reduced chains can 
be systematically compared. This is the representa¬ 
tion described in Fig. [^a). 


2. Convolutions II: Recursively convolving with 
exponential chain 


If we write Eqs. ( |B11[ ) in the form 

(5 -I- ) , 

{s + yTo^s'), 

(C2) 

where Eg = * • • • then we can recursively 

convolute the objects involving S functions onto the 
chain of exponentials from the right. We note that 


2o|zo [tj 

^0|^o 

PNjio (i) 
^N|^o 


= Ejq—i * ((5 -|- ^< 5 ^) * • • • * 
= Ejv— 1 * (<5 + z/i * • • • * 




Eli* [6 ■ 


M'j = 


^E^_i{t)+fl-^]Eii{t) 


(C3) 
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which follows directly from Eq. (A7). From this, each 
of the recursive convolutions introduces a new expo¬ 
nential chain with one step less. For example, 


= En-1 * {5 + + 

^0|2o 


1 - 


= <^ 1 - 


Aat-] 

Xi 

Aat-i 


En-1 H— -—-En-2 

Xi 


*{5 + x^ * ... * (d -h Xfl_^^_^5') 


Xi 

Xn-1 Ajv-2 

Xi X2 


1 - 


Xn-1 


X2 


1 - 


Ajv-1 \ Xn-1 , Xn-1 


Xi 


X2 


En-1 + 

En-3 ^ * (d x^^S') *■■■* {6 + . 


Xl 


1 - 


Xn-2 


X2 


En-2 


(C4) 


Performing all the convolutions leads to the following expressions: 




EN-a{t) I Xn-1: 


9=1 


'y {xji xn-ji) y ^ {xj2+i Ajv-jj) y ^ 

il = l 92=91 


{xjj^_i^_^+N-io-a-l XN-jiv_ig_^) 


93=92 jN-iQ-c = 

jN-ig-a-l 


PNHoit) 1 ^ 


N\io( 

4’N\i 


io-1 

n;ME 


T Vuc / T 

a—1 / a—1 


EN-a{t) j Xn-1: 


E! ivii Xn-Ji) E! (2/92-1-1 Aat-jj) E! ■■■ E! ( 2 / 9 io-Q+*o-a-l AjV-jig-o 

il—1 32—jl j3—j2 jiQ-a = 

jin—a-I 


(C5) 


These expressions can be written as 

N-io 

PolioW = ^0\io E ^N-aP^-ait), 

a=l 

ip 

PN\^o{t) = ^N\^o E GEa^^lV-a(i), (C6) 

a—1 

where P'^N-a constants (independent 

of time). The fixation time distributions are thus 
linear combinations of the distributions En-u- We 
note here the equalities 


states if currently in eigenstate a\ either a —> a -|- 1 
or a ^ N. These paths have transition rate FaXa 
and (1 —Fa)Aa, respectively. The total rate of transi¬ 
tioning out of a is then Aq,, and the waiting time at a 
is an exponential distribution with parameter Aq , no 
matter whether the system transitions to a -I- 1 or to 
N. The quantity Fq denotes the probability that the 
next state of dynamics in eigenspace is a -I- 1, if the 
system is currently in eigenstate a. With probability 
1 — Fq the next state is eigenstate N. Evaluating the 
probability of a tra jecto ry in terms of Fq, and then 
matching with Eq. (C61 gives 


N-ip 

T. = 1. 

a=l 


*0 

E cSv’b = 1 . 

a=l 


(C7) 


We now proceed to express the above linear combi¬ 
nation of exponential chains [Eq. (C 61 ] as the single 
chain shown in Fig. [^b). In the schematic shown in 
this figure the system can transition to two possible 


a 

1 - E G. 

Fq =-El- for a < fV - 1. (C8) 

1 - E G„ 

K — 1 
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We can express all transition rates in Fig. [^b) as 

rp _ f /3 — o “t" 1 < A^, 

1(1-F„)A„, a<lV-l, /3 = W 


distance falls below d = 1/2. We then plot the time 
taken to complete this number of simulation runs of 
the given process. 


3. Evaluation of ‘bottom-line’ arrival time 
distributions 


Appendix D: Relation to equilibration processes 


The final expressions for the (un-normalized ) ar - 
rival time distributions follow directly from Eq. (C6). 
First we note that the convolution of € exponential 
distributions has the form 


Ct—l 0 — 1 

0^a 


Substituting t his expression into Eq. (C6| 
equivalently Eq. (C5l, we arrive at the final ex] 
sions 


or 
expres- 


N-l 


= [UdA E 


C N-1 \ N-1 

E 

£—io J ct—l 


N-io-1 

J~J (^7 

7^1 _ 

N-1 

n (A /3 - Xa) 

13=1 
■ /5/a 

io —1 

n iVl ~ 

7=1 


— Xcct 


-Xc^t 


N-1 

n (A /3 - Aa) 
/ 3=1 
./3^a 


(Cll) 

Evaluating these expressions requires the calculation 
of the eigenvalues Xa,ya, and A^. 


4. Accuracy and efficiency of simulation-based 
approaches 


In the inset of Fig. we show that computing ar¬ 
rival time distributions directly from simulations is 
less efficient than computing them exactly using our 
formalism. To generate the CPU-time curve for the 
simulation method we measure the distance, d, be¬ 
tween the histogram of arrival times at state N, p^, 
against the exact distribution. This distance is de¬ 
fined by 


d 



i=0 



^i) PNidi) 


. + 1 p. 


4’N\io 


df 


(C12) 

where M is the number of histogram bins. Note that 
this is the continuous analog of the distance measure 
between distributions to be discussed in the next sec¬ 
tion. The histogram, pjv, is populated with an in¬ 
creasing amount of independent realizations until the 


1. Dynamics without mutation 


In the system without mutation all realizations 
reach fixation eventually. If the dynamics is started 
from state iq [i.e., Pi(t = 0) = the stationary 
state of the birth-death process, i.e., the terminal dis¬ 
tribution, is of the form 

= (1 ~ + (t)N\iodi,N, (* = 0, . . . , A). 

(Dl) 

The quantity 4)N\io probability that the process 
reaches the absorbing state N. The probability of 
being absorbed at 0 is 1 — (l)N\io- 

Let us now consider the distance of the dis¬ 
tribution P(t) from this distribution $ = (I — 
</ 7 V|io> 0, • • •, 0,/) 7 v|io)- In line with the existing lit¬ 
erature [m [23] we use the distance measure 

1 ^ 

d[P,Q] = -El^*-Q*l’ (D2) 

for two distributions P and Q. We then have 


d[p(t),$] 


PoW - (1 - /’ATlio) 


N-1 

+ E ~ 4>N\io 

i=l 


(D3) 


Probability continuously flows into the absorbing 
states, hence Po(t) < 1 — </Ar|io ^'^d PAr(t) < (l)N\ia 
[Po(f) approaches I — (l>N\io from below with time, 
and similarly for Pn (t)]. We can therefore simplify 
the above expression, and we are left with 


d[P(t),$] = - 


N-1 


I - Po{t) - PN{t) + E P^it) 

= I - Po(t) - PNit). (D4) 


This means that the distance d{t) = d[P(t),$] is 
given by the probability that the system has not yet 
reached fixation in either of the absorbing states by 
time t. This in turn means that l — d{t) = Pr(7x < t) 
is the probability to have reached fixation by time t, 
i.e., it is the cumulative distribution of the uncondi¬ 
tional fixation time tflx- The quantity —d{t) is there¬ 
fore the probability density function of the uncondi¬ 
tional fixation time. 
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As a side remark we note that the mean uncondi¬ 
tional fixation time can be expressed as follows 


(ifix) = 


-d{t) 


Jo 


t dt 


d{t) dt 


d{t) dt. 


(D5) 


Thus the mean unconditional fixation time is the area 
under the curve d{t). 


Together with the normalization condition 
(Sfco = 1) can determine that Pq* and 
P^ must be O (u°), and the remaining probability 
masses are O {u). With this we can write 

Pf = (1 - a)<5,,o + a6,,N + 0{u), (D9) 

for i = 0,..., A^, where a = O (u°). This indicates 
that in the limit of small mutation the distribution is 
peaked at the boundaries. 


c. Intermediate distribution 


2. Dynamics with mutation 

a. Definitions 


We now consider systems with mutation, which oc¬ 
curs with rate u <§; 1. This means that the states 
0 and N are no longer absorbing. More specifi¬ 
cally we consider systems in which bo = O (u) and 
d]^ = 0(u), i.e., escape from the states 0 and N 
occurs with a rate proportional to u. All remaining 
transition rates {bi,di, i = 1,... ,N — 1) are O (m°) , 
and hence are only affected at sub-leading order by 
the introduction of mutation. For u = 0, one recovers 
the case with absorbing states {bo = d^ = 0)- Be¬ 
low we will compare the system with small mutation 
rates with the system without mutation. 

The rates used for the analysis shown in Figs.[^[^ 
and 1^ are given by (see Ref. [53]) 


b^ = 


di = 


(1 - ^^ g[-fA7r(z)] -h 

(1 - ^\ [-ATr{i)] + 


u {N — 
2 N 


• 2 

U 

2N' 


(D6) 


b. Stationary state 


For n > 0 there are no absorbing states and the 
dynamics reaches a stationary distribution, P®*, with 
full support, i.e., Pf* > 0 for alH = 0,1,..., N. This 
distribution can be expressed as EJ 




st 

0 5 




N 


-1 


i+En 


bj-i 


(D7) 


Z=1 j = l 


=1 


In the limit of small mutation rates (0 < uN ^ 1), 
it can be seen that 





0{u) for i = 1 ,..., A^ — 1. 

(D8) 


The system without mutation {u = 0) has termi¬ 
nal distribution €>, as discussed above. In particular 
Po{t) and PN{t), the probability masses concentrated 
in the absorbing states, grow with time and we have 
Po{t) — >■ 1 — 'pN\io PAr(t) —>■ (l>N\io aS t — >■ OO. 

The rates of the system with small mutation differ 
from those of the system without mutation only by 
corrections of O (u). At small mutation rates we ex¬ 
pect the dynamics on a short timescale (t <C to 
be essentially that of the system without mutation, 
the effects of mutation only set in at a longer time. Of 
course the boundary states are no longer absorbing, 
but we argue that the system initially approaches a 
distribution close to $ before reaching its stationary 
distribution P®*. 

This can be seen mathematically as follows: 

Let p(“=o)(t) = [Po^“=°^(t),...,P^"=°^(t)] be the 
probability distribution of the system without muta¬ 
tion. The time evolution is described by the master 
equation 


p(“=o) =M-P(“=°), (DIO) 

where M is an {N -|- 1) x {N + I) transition matrix 
(to be distinguished from the truncated matrix A). 
Let P(“^(t) be the distribution in the system with 
mutation whose evolution is described by 

P(“) = (M-F mQ) • P(“\ (Dll) 

where the matrix Q reflects the difference between the 
systems with and without mutation. The elements of 
both matrices M and Q are O (u^)- Now, let q(t) = 
p(«)(t) _ p(“=o)(t), such that 

q = M-q-hMQ-P(“). (D12) 

We want to calculate how the separation, q, grows 
in time given that both systems (with and without 
mutation) start from the same initial condition [i.e., 
q(0) = 0]. For this purpose it is convenient to work 
in the eigenspace of M. This matrix has two zero 
eigenvalues, no = Di = Oj eigenvectors (vo)i = 
5ifi and (vi)i = These are the absorbing states 
of the system without mutation. 
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Decomposing q(t) = J2a Q{t)a'Va into eigendirec- 
tions Va of M we have 

Qa = + Ugci{t), (D13) 

where we have written Q • ga{t)'Va and 

we note that ga{t) = O (m°)- This can be integrated 
to give 


ga(i)=ne'"“* [ e dr. (D14) 

Jo 

On short time scales {t <C u~^) q(f) = O (u), and 
hence the separation q(t) is also O {u). That is to say 
in the limit u —> 0, both systems (with and without 
mutation) initially evolve along the same trajectory. 
On this time scale both systems approach the distri¬ 
bution 

On a longer time scale [t = O (u“^)], differences 
between the two systems become of O (u°). However 
these differences are concentrated on the states i = 0 
and i = N. Effectively, a redistribution of probability 
mass between the boundary states takes place. The 
distribution of the system with mutation evolves from 
^*|io = (l-'('Af|io)^*.O + 0Ar|io<^bAf toPf = (l-cr)5j,o + 
aSi^Pf, as shown in Fig. 


d. Distance from stationarity and mixing times 

Approximating the stationary distribution of the 
system with small mutation rate as P®* = {l—a)Sifi + 
aSi^N we have 




Pt\t) - (1 - a) 

N-1 


+ E ^ w + 


2=1 


(D15) 


for the distance between the distribution of 

the system with mutation at time t and the station¬ 
ary distribution. While Po(t) and P/v(t) are mono- 
tonically increasing with time in the system without 
mutation, this is not necessarily the case if there is 
mutation. Hence we cannot easily drop the absolute 
values in Eq. (D15l as in the case without mutation. 


We observe, though, that Po’^^(t = 0) = 0 and 
P^\t = 0) = 0 for 0 < *0 < A^. Hence there is an 
initial phase of the dynamics in which Pg“^ (t) < 1 — a 
and P^\t) < a. Let us write t* for the first time at 
which either Po{t*) = 1 — a or P/v(t*) = a (whichever 
happens first). Prior to this time we have 


d[p(“)(t),P'^‘] 


1 

2 


(1 - a) - Po^“)(t) 


N-1 

^=1 


= i_pW(i)_pW(i). (Die) 


This is the same as the distance to the fixation dis¬ 
tribution, €>, in the system without mutation, given 
in Eq. ( |D4[ ). From this we can conclude that 

d[p(“)(t),P"*] = d[p(“)(t),$] for t < t*. (D17) 


This is illustrated in Figs. and where the distri¬ 
butions P®* and $ are represented as single points in 
the (Po,P/v) plane. 

Therefore, for times t < t*, we have the relation 


d[p(“)(<),P"*] = d[p(“)(t),$] 

«d[p(“=°)(f),$] =Pr(tflx>t). 

(D18) 


That is to say, the mixing time is related to the cumu¬ 
lative distribution of fixation times. The first time at 
which d[P(“)(t),P®‘] = e, provided t < t*, is approx¬ 
imately the (1 — e)th percentile of the unconditional 
fixation time distribution. Choosing e = 1/2, we have 
the equivalence with between the mixing time and the 
median fixation time. This equivalence holds for the 
example shown in Figs. and However, it does 
not hold in the example shown in Fig. where the 
stationary state is asymmetric in t he (P p, Py) plane. 
In this case the equivalence in Eq. ( D18[ ) breaks down 
prior to the time at which d[p(“) (t), P®*^] = 1/2. 
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